Análise Estatística Multivariada

1 Análise de Correlação Canônica

O termo correlação se refere ao relacionamento entre variáveis e é medida pelo coeficiente de correlação. A correlação simples mede o grau de associação linear entre duas variáveis (\(r_{xy}\)). A correlação parcial mede o grau de associação linear entre duas variáveis mantendo-se outras variáveis constantes. A ordem da correlação parcial vai depender do número de variáveis que ficam constantes. Com uma variável constante (\(r_{xy.z}\)) Z, a correlação parcial de x e y é dita de ordem um. A correlação múltipla mede o grau de associação entre uma variável e um conjunto de variáveis (\(r_{x_1.x_2 x_3 \dots x_p}\)). A correlação canônica mede o grau de associação entre dois grupos de variáveis, com um grupo denotado normalmente por Y e o outro por X. Estatiscamente, a correlação canônica procura testar se os dois conjuntos de variáveis são independentes.

É uma técnica de Análise Multivariada usada para quantificar o grau de correlação (dependência) entre dois grupos de variáveis. Exemplos:

1. Correlação entre Variáveis Econômicas e Variáveis Sociais
2. Correlação entre Variáveis Tecnológicas e Variáveis Econômicas
3. Correlação entre Variáveis Econômicas e Variáveis Ambientais
4. Correlação entre Variáveis de Habilidade Matemática e Variáveis de Habilidade de Escrita
5. Correlação entre Variáveis de Características do Trabalho e Variáveis de Satisfação no Trabalho

Um exemplo disponível no site da UCLA (https://stats.oarc.ucla.edu/r/dae/canonical-correlation-analysis/) mostra que uma pesquisadora coletou dados sobre três variáveis psicológicas, quatro variáveis acadêmicas (resultados de testes padronizados) e gênero para 600 calouros universitários. Ela está interessada em como o conjunto de variáveis psicológicas se relaciona com as variáveis acadêmicas e de gênero. Em particular, a pesquisadora está interessada em quantas dimensões (variáveis canônicas) são necessárias para entender a associação entre os dois conjuntos de variáveis.

O arquivo de dados (mmreg.csv), possui 600 observações em oito variáveis. As variáveis psicológicas são locus_of_control, auto_conceito e motivação. As variáveis acadêmicas são testes padronizados em leitura (read), escrita (write), matemática (math) e ciências (science). Além disso, a variável female é uma variável binária (zero-um) com 1 para a estudante do sexo feminino.

1.1 Entrada Dados da Correlação Canônica no R

#Pacotes utilizados
library(tidyverse)
library(skimr)
library(ggplot2)
library(GGally)
library(CCA)
library(CCP)
library(psych)

dados <- read.csv("https://stats.idre.ucla.edu/stat/data/mmreg.csv")
colnames(dados) <- c("Control", "Concept", "Motivation", "Read", "Write", "Math", 
    "Science", "Gender")

attach(dados)

#Visualizaçao dos dados
glimpse(dados)

## Rows: 600
## Columns: 8
## $ Control    <dbl> -0.84, -0.38, 0.89, 0.71, -0.64, 1.11, 0.06, -0.91, 0.45, 0…
## $ Concept    <dbl> -0.24, -0.47, 0.59, 0.28, 0.03, 0.90, 0.03, -0.59, 0.03, 0.…
## $ Motivation <dbl> 1.00, 0.67, 0.67, 0.67, 1.00, 0.33, 0.67, 0.67, 1.00, 0.67,…
## $ Read       <dbl> 54.8, 62.7, 60.6, 62.7, 41.6, 62.7, 41.6, 44.2, 62.7, 62.7,…
## $ Write      <dbl> 64.5, 43.7, 56.7, 56.7, 46.3, 64.5, 39.1, 39.1, 51.5, 64.5,…
## $ Math       <dbl> 44.5, 44.7, 70.5, 54.7, 38.4, 61.4, 56.3, 46.3, 54.4, 38.3,…
## $ Science    <dbl> 52.6, 52.6, 58.0, 58.0, 36.3, 58.0, 45.0, 36.3, 49.8, 55.8,…
## $ Gender     <int> 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 1, 0, 1, 1, 1, 1,…

head(dados)

##   Control Concept Motivation Read Write Math Science Gender
## 1   -0.84   -0.24       1.00 54.8  64.5 44.5    52.6      1
## 2   -0.38   -0.47       0.67 62.7  43.7 44.7    52.6      1
## 3    0.89    0.59       0.67 60.6  56.7 70.5    58.0      0
## 4    0.71    0.28       0.67 62.7  56.7 54.7    58.0      0
## 5   -0.64    0.03       1.00 41.6  46.3 38.4    36.3      1
## 6    1.11    0.90       0.33 62.7  64.5 61.4    58.0      1

tail(dados)

##     Control Concept Motivation Read Write Math Science Gender
## 595    0.46    0.03       1.00 52.1  56.7 62.8    47.1      1
## 596    0.94   -0.30       1.00 60.1  67.1 52.4    55.3      1
## 597    0.23    0.03       1.00 65.4  56.7 65.4    58.0      1
## 598    0.46    0.03       1.00 65.4  51.5 61.4    60.7      1
## 599    0.51    0.03       1.00 54.8  54.1 66.4    41.7      1
## 600    0.25    0.03       0.67 49.5  51.5 55.5    44.4      1

#Estatistica descritiva dos dados
summary(dados)

##     Control            Concept            Motivation          Read     
##  Min.   :-2.23000   Min.   :-2.620000   Min.   :0.0000   Min.   :28.3  
##  1st Qu.:-0.37250   1st Qu.:-0.300000   1st Qu.:0.3300   1st Qu.:44.2  
##  Median : 0.21000   Median : 0.030000   Median :0.6700   Median :52.1  
##  Mean   : 0.09653   Mean   : 0.004917   Mean   :0.6608   Mean   :51.9  
##  3rd Qu.: 0.51000   3rd Qu.: 0.440000   3rd Qu.:1.0000   3rd Qu.:60.1  
##  Max.   : 1.36000   Max.   : 1.190000   Max.   :1.0000   Max.   :76.0  
##      Write            Math          Science          Gender     
##  Min.   :25.50   Min.   :31.80   Min.   :26.00   Min.   :0.000  
##  1st Qu.:44.30   1st Qu.:44.50   1st Qu.:44.40   1st Qu.:0.000  
##  Median :54.10   Median :51.30   Median :52.60   Median :1.000  
##  Mean   :52.38   Mean   :51.85   Mean   :51.76   Mean   :0.545  
##  3rd Qu.:59.90   3rd Qu.:58.38   3rd Qu.:58.65   3rd Qu.:1.000  
##  Max.   :67.10   Max.   :75.50   Max.   :74.20   Max.   :1.000

skim(dados)

Data summary

Name	dados
Number of rows	600
Number of columns	8
_______________________
Column type frequency:
numeric	8
________________________
Group variables	None

Variable type: numeric

skim_variable	n_missing	complete_rate	mean	sd	p0	p25	p50	p75	p100	hist
Control	0	1	0.10	0.67	-2.23	-0.37	0.21	0.51	1.36	▁▂▃▇▃
Concept	0	1	0.00	0.71	-2.62	-0.30	0.03	0.44	1.19	▁▁▃▇▃
Motivation	0	1	0.66	0.34	0.00	0.33	0.67	1.00	1.00	▂▃▁▆▇
Read	0	1	51.90	10.10	28.30	44.20	52.10	60.10	76.00	▂▇▇▆▂
Write	0	1	52.38	9.73	25.50	44.30	54.10	59.90	67.10	▁▃▅▆▇
Math	0	1	51.85	9.41	31.80	44.50	51.30	58.37	75.50	▃▇▇▅▂
Science	0	1	51.76	9.71	26.00	44.40	52.60	58.65	74.20	▁▅▆▇▂
Gender	0	1	0.54	0.50	0.00	0.00	1.00	1.00	1.00	▇▁▁▁▇

Abaixo está uma lista de alguns métodos de análise que podem ser usadas com estes dados:

1. Análise de correlação canônica;

2. Regressões simples separadas – analisar esses dados usando análises de regressão simples separadas para cada variável em um conjunto. As regressões não produzirão resultados multivariados e não relatam informações relativas à dimensionalidade;

3. A regressão múltipla multivariada.

Os dados para análise de correlações canônicas são formados por dois conjuntos de variáveis, formando dois vetores de variáveis X e Y, que podem ser organizados como demonstrado abaixo.

OBS			X				Y
	X1	X2	...	XP	Y1	Y2	...	YP
1
2
3
4
...
N

com estes dados pode-se ter uma matriz de correlaçoes dada por

\[ \mathbf{R}=\left[\begin{array}{ccc} R_{11} & \vdots & R_{12} \\ \dots & \vdots & \dots \\ R_{21} & \vdots & R_{22} \\ \end{array} \right] \]

que pode ser particionada em 4 partes onde \(R_{11}\) é uma matriz de correlações do primeiro grupo; \(R_{22}\) é uma matriz de correlações do segundo grupo; e \(R_{12}\) e \(R_{21}\) são correlações entre as variáveis de um grupo com as do outro grupo, sendo que uma é a transposta da outra.

1.2 Verificação das Correlações no R

psych <- dados[, 1:3]
acad <- dados[, 4:8]

ggpairs(psych)

ggpairs(acad)

# Correlações entre os dois conjuntos de variáveis
matcor(psych, acad)

## $Xcor
##              Control   Concept Motivation
## Control    1.0000000 0.1711878  0.2451323
## Concept    0.1711878 1.0000000  0.2885707
## Motivation 0.2451323 0.2885707  1.0000000
## 
## $Ycor
##                Read     Write       Math    Science      Gender
## Read     1.00000000 0.6285909  0.6792757  0.6906929 -0.04174278
## Write    0.62859089 1.0000000  0.6326664  0.5691498  0.24433183
## Math     0.67927568 0.6326664  1.0000000  0.6495261 -0.04821830
## Science  0.69069291 0.5691498  0.6495261  1.0000000 -0.13818587
## Gender  -0.04174278 0.2443318 -0.0482183 -0.1381859  1.00000000
## 
## $XYcor
##              Control     Concept Motivation        Read      Write       Math
## Control    1.0000000  0.17118778 0.24513227  0.37356505 0.35887684  0.3372690
## Concept    0.1711878  1.00000000 0.28857075  0.06065584 0.01944856  0.0535977
## Motivation 0.2451323  0.28857075 1.00000000  0.21060992 0.25424818  0.1950135
## Read       0.3735650  0.06065584 0.21060992  1.00000000 0.62859089  0.6792757
## Write      0.3588768  0.01944856 0.25424818  0.62859089 1.00000000  0.6326664
## Math       0.3372690  0.05359770 0.19501347  0.67927568 0.63266640  1.0000000
## Science    0.3246269  0.06982633 0.11566948  0.69069291 0.56914983  0.6495261
## Gender     0.1134108 -0.12595132 0.09810277 -0.04174278 0.24433183 -0.0482183
##                Science      Gender
## Control     0.32462694  0.11341075
## Concept     0.06982633 -0.12595132
## Motivation  0.11566948  0.09810277
## Read        0.69069291 -0.04174278
## Write       0.56914983  0.24433183
## Math        0.64952612 -0.04821830
## Science     1.00000000 -0.13818587
## Gender     -0.13818587  1.00000000

A matriz \(R_{12}\) mostra as correlações entre as variáveis dos dois grupos. Com p variáveis no vetor X e q variáveis no vetor Y tem-se pq correlações. A análise do relacionamento entre os dois grupos de variáveis por meio das pq correlações é inviável. O objetivo da Análise de Correlação Canônica é resumir o padrão de associação entre os X’s e os Y’s em termos de poucas correlações.

A idéia básica da técnica é:

1. Calcula-se, inicialmente, 2 combinações lineares, uma de cada conjunto de variáveis, de forma que a correlação entre elas seja máxima. As combinações lineares são denominadas variáveis canônicas e os pares são denominados pares canônicos. O primeiro par deve apresentar a maior correlação entre as variáveis canônicas;

2. Em seguida, calcula-se 2 outras variáveis canônicas, ou seja, outro par canônico, com a condição de serem ortogonais às primeiras e terem o máximo coeficiente de correlação. O segundo par deve ter a segunda maior correlação entre as variáveis canônicas.

3. Assim, sucessivamente.

Em geral, o número de dimensões canônicas (pares canônicos) é igual ao número de variáveis no conjunto menor \(min(p,q)\). No entanto, o número de dimensões significativas pode ser ainda menor. O coeficiente de correlação canônica é o coeficiente de correlação de Pearson, em valor absoluto, entre as variáveis canônicas de cada par canônico. Existem, então, \(k=min(p,q)\) pares canônicos e k coeficientes de correlação canônica sendo interpretados somente aqueles estatísticamente significativos.

Sejam \(U_1\) e \(V_1\) o primeiro par canônico, dado por

\[\begin{cases} U_1=a_{11}X_1+a_{21}X_2+a_{31}X_3+ \dots +a_{p1}X_p \\ V_1=b_{11}Y_1+b_{21}Y_2+b_{31}Y_3+ \dots +b_{q1}Y_q \\ \end{cases}\]

o problema passa a ser encontrar os vetores \(a_1\) e \(b_1\) que maximizem a correlação entre \(U_1\) e \(V_1\).

Sejam \(U_2\) e \(V_2\) o segundo par canônico, dado por

\[\begin{cases} U_2=a_{12}X_1+a_{22}X_2+a_{32}X_3+ \dots +a_{p2}X_p \\ V_2=b_{12}Y_1+b_{22}Y_2+b_{32}Y_3+ \dots +b_{q2}Y_q \\ \end{cases}\]

o problema passa a ser agora encontrar os vetores \(a_2\) e \(b_2\) que maximizem a correlação entre \(U_2\) e \(V_2\) e mantenha a condição de ortogonalidade entre \(U_1\) e \(U_2\) e entre \(V_1\) e \(V_2\).

1.3 Solução Geral da Correlação Canônica

Considerando as variáveis padronizadas, que é o mesmo que trabalhar com a matriz de correlações, tem-se que combinações lineares das variáveis X são dadas por \(U=a'X\), com \(var(U)=var(a'X)=a'R_{11}a\) e combinações lineares de Y são dadas por \(V=b'Y\), com \(var(V)=var(b'Y)=b'R_{22}b\) e \(cov(U,Y)=cov(a'X,b'Y)=a'R_{12}b\).

Sendo

\[ r_{UV}=\frac{COV(U,V)}{\sqrt{var(U)}\sqrt{var(V)}} \]

pode-se definir \(r_{UV}^2\) como

\[ r_{UV}^2=\frac{[COV(U,V)]^2}{var(U) var(V)}=\frac{(a'R_{12}b)^2}{(a'R_{11}a)(b'R_{22}b)} \]

Assim, o problema pode ser definido como

\[ maximizar \quad (a'R_{12}b) \]

sujeito a

\[ a'R_{11}a=1 \] \[ b'R_{22}b=1 \] \[ a'R_{11}a^*=0 \quad a \neq a^* \] \[ b'R_{22}b^*=0 \quad b \neq b^* \]

Resumidamente, o problema de maximização condicionado encontra os valores de \(a_k\) e \(b_k\) como soluções do seguinte sistema de equações

\[\begin{cases} (R_{12}R_{22}^{-1}R_{21}-\lambda_k R_{11})a_k=0 \\ (R_{21}R_{11}^{-1}R_{12}-\lambda_k R_{22})b_k=0 \\ \end{cases}\]

em que \(\lambda_k\) são raízes características e \(a_k\) e \(b_k\) são vetores característicos. Assim, os \(\lambda_k\) satisfazem as seguintes equações características

\[\begin{cases} |R_{12}R_{22}^{-1}R_{21}-\lambda_k R_{11}|=0 \\ |R_{21}R_{11}^{-1}R_{12}-\lambda_k R_{22}|=0 \\ \end{cases}\]

ou

\[\begin{cases} |R_{11}^{-1}R_{12}R_{22}^{-1}R_{21}-\lambda_k I|=0 \\ |R_{22}^{-1}R_{21}R_{11}^{-1}R_{12}-\lambda_k I|=0 \\ \end{cases}\]

em que, \(\lambda_k\) é a k-ésima maior raiz característica da matriz \(R_{11}^{-1}R_{12}R_{22}^{-1}R_{21}\) ou da matriz \[R_{22}^{-1}R_{21}R_{11}^{-1}R_{12}\].

O vetor “a” é o vetor característico da matriz \(R_{11}^{-1}R_{12}R_{22}^{-1}R_{21}\)

Dado “a”, o vetor “b” pode ser obtido pela relação \(b=\frac{1}{\sqrt{\lambda}}R_{22}^{-1}R_{21}a\).

O vetor “b” pode ser obtido pelas raízes características da matriz \(R_{22}^{-1}R_{21}R_{11}^{-1}R_{12}\).

De forma semelhante, dado “b”, o vetor “a” pode ser calculado por \(a=\frac{1}{\sqrt{\lambda}}R_{11}^{-1}R_{12}b\).

As soluções para as raízes de \(R_{11}^{-1}R_{12}R_{22}^{-1}R_{21}\) ou de \[R_{22}^{-1}R_{21}R_{11}^{-1}R_{12}\] são idênticas.

Tem-se, assim, a solução para os coeficientes das variáveis canônicas que são os elementos dos vetores “a” e “b”. Observe que:

1. As raízes podem ser calculadas com base em equações características diferentes

\[\begin{cases} |R_{11}^{-1}R_{12}R_{22}^{-1}R_{21}-\lambda_k I|=0 \\ |R_{22}^{-1}R_{21}R_{11}^{-1}R_{12}-\lambda_k I|=0 \\ \end{cases}\]

2. Para cada raiz existe um par de vetores cujos elementos são os coeficientes das variáveis canônicas. Estes coeficientes são chamados pesos canônicos e podem ser calculados na forma bruta ou padronizada.

3. O coeficiente de correlação canônica é o coeficiente de correlação, em valor absoluto, entre \(U_k\) e \(V_k\). Ele é igual a raiz quadrada da raiz característica, ou seja, \(r_1=\sqrt{\lambda_1}\), \(r_2=\sqrt{\lambda_2}\), etc.

1.4 Análise da Correlação Canônica no R

# Calcular e mostrar as correlaçoes canônicas
cc1 <- cc(psych, acad)
cc1$cor #mostra as correlaçoes

## [1] 0.4640861 0.1675092 0.1039911

cc1[3:4] #mostra os coeficientes canônicos brutos

## $xcoef
##                  [,1]       [,2]       [,3]
## Control    -1.2538339 -0.6214776 -0.6616896
## Concept     0.3513499 -1.1876866  0.8267210
## Motivation -1.2624204  2.0272641  2.0002283
## 
## $ycoef
##                 [,1]         [,2]         [,3]
## Read    -0.044620600 -0.004910024  0.021380576
## Write   -0.035877112  0.042071478  0.091307329
## Math    -0.023417185  0.004229478  0.009398182
## Science -0.005025152 -0.085162184 -0.109835014
## Gender  -0.632119234  1.084642326 -1.794647036

A dimensão 1 teve uma correlação canônica de 0,46 entre os conjuntos de variáveis, enquanto para a dimensão 2 a correlação canônica foi bem menor em 0,17. A correlação canônica da terceira dimensão foi 0,10.

Os coeficientes canônicos brutos são interpretados da mesma forma que os coeficientes da análise de regressão, ou seja, para a variável read, um aumento de uma unidade leva a uma diminuição de 0,0446 na primeira variável canônica do segundo conjunto de dados, tudo o mais constante. Outro exemplo: ser mulher leva a uma diminuição de 0,6321 na dimensão 1 para o conjunto de dados “acadêmico” com tudo mais constante.

# Calcular as cargas canônicas
cc2 <- comput(psych, acad, cc1)

# Mostrar as cargas canônicas
cc2[3:6]

## $corr.X.xscores
##                   [,1]       [,2]       [,3]
## Control    -0.90404631 -0.3896883 -0.1756227
## Concept    -0.02084327 -0.7087386  0.7051632
## Motivation -0.56715106  0.3508882  0.7451289
## 
## $corr.Y.xscores
##               [,1]        [,2]        [,3]
## Read    -0.3900402 -0.06010654  0.01407661
## Write   -0.4067914  0.01086075  0.02647207
## Math    -0.3545378 -0.04990916  0.01536585
## Science -0.3055607 -0.11336980 -0.02395489
## Gender  -0.1689796  0.12645737 -0.05650916
## 
## $corr.X.yscores
##                    [,1]        [,2]        [,3]
## Control    -0.419555307 -0.06527635 -0.01826320
## Concept    -0.009673069 -0.11872021  0.07333073
## Motivation -0.263206910  0.05877699  0.07748681
## 
## $corr.Y.yscores
##               [,1]        [,2]       [,3]
## Read    -0.8404480 -0.35882541  0.1353635
## Write   -0.8765429  0.06483674  0.2545608
## Math    -0.7639483 -0.29794884  0.1477611
## Science -0.6584139 -0.67679761 -0.2303551
## Gender  -0.3641127  0.75492811 -0.5434036

As correlações acima são entre variáveis observadas e variáveis canônicas que são conhecidas como cargas canônicas. Essas variáveis canônicas são na verdade um tipo de variável latente, ou seja, não observadas. À semelhança da interpretação dos fatores na Análise Fatorial, as variáveis canônicas podem ser identificadas em termos das variáveis com que elas mais se relacionam. Isto pode ser feito por meio dos coeficientes canônicos ou pesos canônicos ou, preferencialmente, através das correlações entre a variável canônica e as variáveis originais. Estas correlações são denominadas cargas canônicas (canonical loadings) ou correlações estruturais. Por exemplo, se \(U_1\) e \(X_1\) são positiva e altamente correlacionadas, \(U_1\) pode ser considerada um reflexo de \(X_1\). Da mesma forma, se \(V_1\) for negativa e altamente correlacionada com \(Y_1\) e positiva e altamente correlacionada com \(Y_3\), então \(V_1\) reflete \(Y_1\) em sentido contrário e \(Y_3\) no mesmo sentido.

São calculados coeficientes de correlação entre as variáveis canônicas e as variáveis originais padronizadas de forma que se tenha os coeficientes de correlação entre as variáveis canônicas U e as variáveis do grupo X; entre as variáveis canônicas U e as variáveis do grupo Y; entre as variáveis canônicas V e as variáveis do grupo Y e entre as variáveis canônicas V e as variáveis do grupo X.

Quando as variáveis do modelo possuem desvios padrão muito diferentes, os coeficientes padronizados permitem comparações mais fáceis entre as variáveis, calculados como demonstrado abaixo.

# coeficientes canonicos padronizados para  psych (desvio padrao da diagonal da matriz de var-cov de psych)
s1 <- diag(sqrt(diag(cov(psych))))
s1 %*% cc1$xcoef

##            [,1]       [,2]       [,3]
## [1,] -0.8404196 -0.4165639 -0.4435172
## [2,]  0.2478818 -0.8379278  0.5832620
## [3,] -0.4326685  0.6948029  0.6855370

# coeficientes canonicos padronizados para  acad (desvio padrao da diagonal da matriz de var-cov de acad)
s2 <- diag(sqrt(diag(cov(acad))))
s2 %*% cc1$ycoef

##             [,1]        [,2]        [,3]
## [1,] -0.45080116 -0.04960589  0.21600760
## [2,] -0.34895712  0.40920634  0.88809662
## [3,] -0.22046662  0.03981942  0.08848141
## [4,] -0.04877502 -0.82659938 -1.06607828
## [5,] -0.31503962  0.54057096 -0.89442764

Os coeficientes canônicos padronizados são interpretados de maneira semelhante à interpretação dos coeficientes de analise de regressão com variáveis padronizadas. Por exemplo, considere a variável read, um aumento de um desvio padrão em read leva a uma diminuição de 0,45 desvio padrão na pontuação na primeira variável canônica para o segundo conjunto de dados, tudo mais constante.

Para as variáveis psicológicas, a primeira dimensão canônica é mais fortemente influenciada pelo locus de controle (Control) (-0,84) e, para a segunda dimensão, autoconceito (Concept) (-0,84) e motivação (Motivation) (0,69). Para as variáveis acadêmicas mais gênero, a primeira dimensão foi composta por leitura (reading) (-0,45), escrita (writing) (-0,35) e gênero (gender) (-0,32). Para a segunda dimensão, escrita (writing) (0,41), ciência (science) (-0,83) e gênero (gender) (0,54) foram as variáveis dominantes.

Calculados os vetores de coeficientes a e b e os coeficientes de correlação canônica os próximos passos consistem em realizar e interpretar os testes de hipóteses sobre o número de dimensões canônicas. As fórmulas para cálculo podem ser encontradas nas referencias do pacote CCP (https://cran.r-project.org/web/packages/CCP/CCP.pdf)

# testes de dimensões canonicas
rho <- cc1$cor

## Define o numero de observações, numero de variaveis do primeiro grupo e de variáveis do segundo conjunto de variáveis.

n <- dim(psych)[1]
p <- length(psych)
q <- length(acad)

## Calcula os p-values de diferentes testes estatísticos:
p.asym(rho, n, p, q, tstat = "Wilks")

## Wilks' Lambda, using F-approximation (Rao's F):
##               stat    approx df1      df2     p.value
## 1 to 3:  0.7543611 11.715733  15 1634.653 0.000000000
## 2 to 3:  0.9614300  2.944459   8 1186.000 0.002905057
## 3 to 3:  0.9891858  2.164612   3  594.000 0.091092180

p.asym(rho, n, p, q, tstat = "Hotelling")

##  Hotelling-Lawley Trace, using F-approximation:
##                stat    approx df1  df2     p.value
## 1 to 3:  0.31429738 12.376333  15 1772 0.000000000
## 2 to 3:  0.03980175  2.948647   8 1778 0.002806614
## 3 to 3:  0.01093238  2.167041   3 1784 0.090013176

p.asym(rho, n, p, q, tstat = "Pillai")

##  Pillai-Bartlett Trace, using F-approximation:
##                stat    approx df1  df2     p.value
## 1 to 3:  0.25424936 11.000571  15 1782 0.000000000
## 2 to 3:  0.03887348  2.934093   8 1788 0.002932565
## 3 to 3:  0.01081416  2.163421   3 1794 0.090440474

Conforme mostrado na tabela acima, o primeiro teste das dimensões canônicas testa se todas as três dimensões são significativas (elas são, F = 11,72). O próximo teste verifica se as dimensões 2 e 3 combinadas são significativas (elas são, F = 2,94). Por fim, o último teste analisa se a dimensão 3, por si só, é significativa (não é). Portanto, as dimensões 1 e 2 devem ser significativas, enquanto a dimensão três não é.

Um outro exemplo é de dados de uma amostra de 50 funcionários de uma firma que procura aperfeiçoar testes que podem revelar o potencial para vendas. Foram avaliadas 3 medidas de performance: \(X_1\)= aumento de vendas; \(X_2\)= rendimento das vendas, e \(X_3\)= vendas para novos clientes. Essas medidas foram transformadas para uma escala em que 100 indica performance média. Cada indivíduo foi submetido a 4 testes com o objetivo de medir criatividade, raciocínio mecânico, raciocínio abstrato e habilidade matemática. As notas destes testes formam as variáveis: \(Y_1\)= nota no teste de criatividade; \(Y_2\)= nota no teste de raciocínio mecânico; \(Y_3\)= nota no teste de raciocínio abstrato, e \(Y_4\)= nota no teste de habilidade matemática (Johnson e Wichern, 2007, p. 536). Este exemplo é semelhante ao de Mingoti (2005, p. 146).

# Pacotes
library(tidyverse)

# Direcionado o R para o Diretorio a ser trabalhado
setwd('/Users/jricardofl/Dropbox/tempecon/multivariada')

# Entrada dos dados
dados <- read.csv2("vendedor.csv", header=TRUE, sep=";", dec=".")
dados <- dados[,-1] #retira a primeira coluna

#attach(dados)

#Visualizaçao dos dados
glimpse(dados)

## Rows: 50
## Columns: 7
## $ X1 <dbl> 93.0, 88.8, 95.0, 101.3, 102.0, 95.8, 95.5, 110.8, 102.8, 106.8, 10…
## $ X2 <dbl> 96.0, 91.8, 100.3, 103.8, 107.8, 97.5, 99.5, 122.0, 108.3, 120.5, 1…
## $ X3 <dbl> 97.8, 96.8, 99.0, 106.8, 103.0, 99.3, 99.0, 115.3, 103.8, 102.0, 10…
## $ Y1 <int> 9, 7, 8, 13, 10, 10, 9, 18, 10, 14, 12, 10, 16, 8, 13, 7, 11, 11, 5…
## $ Y2 <int> 12, 10, 12, 14, 15, 14, 12, 20, 17, 18, 17, 18, 17, 10, 10, 9, 12, …
## $ Y3 <int> 9, 10, 9, 12, 12, 11, 9, 15, 13, 11, 12, 8, 11, 11, 8, 5, 11, 11, 1…
## $ Y4 <int> 20, 15, 26, 29, 32, 21, 25, 51, 31, 39, 32, 31, 34, 34, 34, 16, 32,…

head(dados)

##      X1    X2    X3 Y1 Y2 Y3 Y4
## 1  93.0  96.0  97.8  9 12  9 20
## 2  88.8  91.8  96.8  7 10 10 15
## 3  95.0 100.3  99.0  8 12  9 26
## 4 101.3 103.8 106.8 13 14 12 29
## 5 102.0 107.8 103.0 10 15 12 32
## 6  95.8  97.5  99.3 10 14 11 21

tail(dados)

##       X1    X2    X3 Y1 Y2 Y3 Y4
## 45  94.8 101.8  99.8  7 16 11 24
## 46 103.5 112.0 110.8 18 13 12 37
## 47  89.5  96.0  97.3  7 15 11 14
## 48  84.3  89.8  94.3  8  8  8  9
## 49 104.3 109.5 106.5 14 12 12 36
## 50 106.0 118.5 105.0 12 16 11 39

#Estatistica descritiva dos dados
summary(dados)

##        X1               X2              X3               Y1       
##  Min.   : 81.50   Min.   : 87.3   Min.   : 94.30   Min.   : 1.00  
##  1st Qu.: 93.55   1st Qu.: 99.5   1st Qu.: 99.08   1st Qu.: 8.25  
##  Median :100.65   Median :106.2   Median :103.15   Median :10.00  
##  Mean   : 98.84   Mean   :106.6   Mean   :102.81   Mean   :11.22  
##  3rd Qu.:105.05   3rd Qu.:114.8   3rd Qu.:106.45   3rd Qu.:14.00  
##  Max.   :110.80   Max.   :122.3   Max.   :115.30   Max.   :18.00  
##        Y2              Y3              Y4       
##  Min.   : 5.00   Min.   : 5.00   Min.   : 9.00  
##  1st Qu.:12.00   1st Qu.: 9.00   1st Qu.:21.50  
##  Median :15.00   Median :11.00   Median :31.50  
##  Mean   :14.18   Mean   :10.56   Mean   :29.76  
##  3rd Qu.:17.00   3rd Qu.:12.00   3rd Qu.:37.00  
##  Max.   :20.00   Max.   :15.00   Max.   :51.00

skim(dados)

Data summary

Name	dados
Number of rows	50
Number of columns	7
_______________________
Column type frequency:
numeric	7
________________________
Group variables	None

Variable type: numeric

skim_variable	n_missing	complete_rate	mean	sd	p0	p25	p50	p75	p100	hist
X1	0	1	98.84	7.34	81.5	93.55	100.65	105.05	110.8	▂▅▃▇▆
X2	0	1	106.62	10.12	87.3	99.50	106.25	114.75	122.3	▅▆▇▇▇
X3	0	1	102.81	4.71	94.3	99.08	103.15	106.45	115.3	▃▅▇▃▁
Y1	0	1	11.22	3.95	1.0	8.25	10.00	14.00	18.0	▁▃▇▅▃
Y2	0	1	14.18	3.38	5.0	12.00	15.00	17.00	20.0	▁▃▆▇▃
Y3	0	1	10.56	2.14	5.0	9.00	11.00	12.00	15.0	▂▃▇▆▁
Y4	0	1	29.76	10.54	9.0	21.50	31.50	37.00	51.0	▅▃▇▆▂

#Correlações entre os dados
perform <- dados %>% 
  dplyr::select(X1, X2, X3) %>%
  scale() %>% as_tibble

testes <- dados %>% 
  dplyr::select(Y1, Y2, Y3, Y4) %>%
  scale() %>% as_tibble

ggpairs(perform)

ggpairs(testes)

matcor(perform, testes)

## $Xcor
##           X1        X2        X3
## X1 1.0000000 0.9260758 0.8840023
## X2 0.9260758 1.0000000 0.8425232
## X3 0.8840023 0.8425232 1.0000000
## 
## $Ycor
##           Y1        Y2        Y3        Y4
## Y1 1.0000000 0.5907360 0.1469074 0.4126395
## Y2 0.5907360 1.0000000 0.3859502 0.5745533
## Y3 0.1469074 0.3859502 1.0000000 0.5663721
## Y4 0.4126395 0.5745533 0.5663721 1.0000000
## 
## $XYcor
##           X1        X2        X3        Y1        Y2        Y3        Y4
## X1 1.0000000 0.9260758 0.8840023 0.5720363 0.7080738 0.6744073 0.9273116
## X2 0.9260758 1.0000000 0.8425232 0.5415080 0.7459097 0.4653880 0.9442960
## X3 0.8840023 0.8425232 1.0000000 0.7003630 0.6374712 0.6410886 0.8525682
## Y1 0.5720363 0.5415080 0.7003630 1.0000000 0.5907360 0.1469074 0.4126395
## Y2 0.7080738 0.7459097 0.6374712 0.5907360 1.0000000 0.3859502 0.5745533
## Y3 0.6744073 0.4653880 0.6410886 0.1469074 0.3859502 1.0000000 0.5663721
## Y4 0.9273116 0.9442960 0.8525682 0.4126395 0.5745533 0.5663721 1.0000000

2 Análise de Agrupamentos (Cluster)

Segundo Mingoti (2005, p.155), “a Análise de Agrupamentos, também conhecida como Análise de Conglomerados, Classificação ou Cluster, tem como objetivo dividir os elementos da amostra, ou população, em grupos de forma que os elementos pertencentes a um mesmo grupo sejam similares entre si com respeito às variáveis (características) que neles foram medidas, e os elementos em grupos diferentes sejam heterogêneos em relação a estas mesmas características”.

Assim, a idéia básica é maximizar a homogeneidade dentro dos grupos ao mesmo tempo em que se maximiza a heterogeneidade entre os grupos. Para formar grupos de objetos tem-se que medir e comparar a semelhança entre eles. A medida de semelhança fica dificultada quando se considera várias variáveis (características).

No processo de agrupamento, primeiro o conjunto de dados será dividido em grupos com elementos parecidos entre si e diferentes dos outros grupos. Depois os grupos encontrados serão analisados até que efetivamente se encontre a existência de padrões entre eles.

Etapas da Análise de Agrupamento - Clustering

No caso univariado é possível que apenas uma análise visual ou gráfica dos dados seja suficiente para decidir sobre os agrupamentos de observaçoes. Contudo, isto não é muito mais difícil quando se tem um conjunto maior de variáveis. Para realizar as técnicas de cluster é necessário definir o conceito de distância.

2.1 Etapas da Análise de Cluster

1. Escolher um critério de parecença (similaridade ou dissimilaridade);

2. Formação dos grupos (escolher o algoritmo de agrupamento);

3. Definição de número de grupos: pode ser a posteriori, como resultado da Análise ou a priori, dado o conhecimento ou conveniência da Análise;

4. Validação do agrupamento: critério do pesquisador;

5. Interpretação e Análise: pode-se fazer estatística descritiva, diferença de médias, etc.

2.2 Medidas de Parecença

Valor numérico que quantifica o grau de semelhança entre um par de objetos. Pode ser de dois tipos: similaridade ou dissimilaridade. Para variáveis quantitativas, quando se buscas agrupar observações as distâncias são as mais usadas. São medidas de dissimilaridade entre objetos. O coeficiente de distância assume valor máximo para objetos totalmente diferentes e valor zero para dois objetos idênticos considerando todas as variáveis.

Em algumas aplicações deseja-se agrupar variáveis ao invés de observações. Neste caso, o coeficiente de correlação de Pearson é mais comumente utilizado. No caso de variáveis binárias, os dados são novamente rearranjados na forma de uma tabela de contigência, com as variáveis formando as categorias.

Os coeficientes de distância mais comuns têm as seguintes propriedades:

1. Mínimo Zero: Se \(A=B\), então \(D(A,B)=0\);

2. Positividade: Se \(A \neq B\), então \(D(A,B)>0\);

3. Simetria: \(D(A,B) = D(B,A)\)

4. Desigualdade triangular: \(D(A,B)+D(B,C) \geq D(A,C)\)

Em relação às distâncias sabe-se pelo Teorema de Pitágoras que \(c^2=a^2+b^2\), assim, \(c=\sqrt{a^2+b^2}\), que é a Distância Euclidiana entre dois pontos (A e B). Em termos de variáveis, tem-se:

\[ D_{AB}=\sqrt{(X_{2A}-X_{2B})^2+(X_{1B}-X_{1A})^2} \] matricialmente, a Distância Euclidiana é definida por \(D_{AB}=\sqrt{(X_a-X_b)'(X_a-X_b)}\) e possui as seguintes propriedades:

1. Base geométrica bem definida;

2. Invariante com relação à transformação de origem;

3. Invariante com relação à transformação ortogonal;

4. Não invariante com relação à transformação de escala;

5. Não invariante com relação à transformação não ortogonal.

A Distância Euclidiana Quadrática é dada por

\[ D_{AB}^2=\sum_{j=1}^p(x_{ja}-x_{jb})^2 \]

A Distância Euclidiana Ponderada é definida por

\[ D_{AB}=\sqrt{(X_a-X_b)'A(X_a-X_b)} \]

com A sendo uma matriz de ponderação (pesos diferentes de acordo com a variância das variáveis). Se \(A=I\), \(D_{AB}\) é a Distância Euclidiana. Se \(A=S^{-1}\), ou seja, inverso da matriz var-cov, tem-se a Distância de Mahalanobis.

Se \(A=diag(S_i^2)^{-1}\), em que \(S_i^2\) é a variância amostral da i-ésima variável aleatória, leva-se em consideração na ponderação apenas as diferenças de variâncias das variáveis.

Se \(A=diag(1/p)\) tem-se a Distância Euclidiana Média;

A Distância de Minkowsky é dada por

\[ D_{AB}=\Big [\sum_i^p w_i | X_{ia} - X_{ib}| ^\lambda \Big]^{\frac{1}{\lambda}} \]

onde \(w_i\) são pesos de ponderação para as variáveis. Se \(\lambda=2\), tem-se a Distância Euclidiana. Esta distância é menos afetada pela presença de valores discrepantes numa amostra do que a distância Euclidiana. Se \(\lambda=1\), tem-se a Distância city-block também conhecida como Manhattan.

Para observações representadas através de variáveis qualitativas, é necessária a criação de variáveis binárias, as quais assumem o valor 1 se uma característica de interesse está presente, e 0, caso contrário. Dessa forma, para um par de observações (i; k) medidos através de p variáveis binárias, considere a seguinte Tabela de Contigência:

		Observação k
		1	0	Total
Observação i	1	a	b	a+b
	0	c	d	c+d
	Total	a+c	b+d	p=a+b+c+d

com base nesta tabela é possível listar algumas medidas de similaridade e a idéia básica de cada uma delas.

\(\frac{a+d}{p}\)	Pesos iguais para pares 1-1 e 0-0
\(\frac{2(a+d)}{2(a+d)+b+c}\)	Peso em dobro para pares 1-1 e 0-0
\(\frac{a+d}{a+d+2(b+c)}\)	Peso em dobro para pares descasados
\(\frac{a}{p}\)	Sem pares 0-0 no numerador
\(\frac{a}{a+b+c}\)	Sem pares 0-0 no numerador e no denominador
\(\frac{2a}{2a+b+c}\)	Sem pares 0 - 0 no numerador e no denominador.
	Peso em dobro para pares 1 - 1.
\(\frac{a}{a+2(b+c)}\)	Sem pares 0 - 0 no numerador e no denominador.
	Peso em dobro para pares descasados.
\(\frac{a}{b+c}\)	Razão entre pares 1 - 1 e pares descasados.

Organiza-se a tabela de contingência e calcula-se o coeficiente de similaridade para cada par de observações. Os coeficientes são dispostos em uma matriz simétrica n x n denominada matriz de similaridade.

2.3 Técnicas para construção de Conglomerados (Clusters)

Para Mingoti (2005, p.164), “as técnicas de conglomerados ou clusters são frequentemente classificadas em dois tipos: técnicas hierárquicas e não hierárquicas, sendo que as hierárquicas são classificadas em aglomerativas e divisivas”.

Os Métodos hierárquicos aglomerativos partem do princípio de que no início do processo de agrupamento têm-se n conglomerados (cada elemento é um conglomerado). Em cada passo do algoritmo, os elementos vão sendo agrupados, formando novos conglomerados até o ponto em que todos os elementos estão num único grupo.

Em termos de variabilidade, no estágio inicial é a menor possível e no final, a máxima possível, pois todos os elementos estão agrupados em apenas 1 cluster. Se dois elementos aparecem juntos num mesmo grupo em algum estçgio do processo de agrupamento, permanecerão juntos até o final, ou seja, não podem mais ser separados.

Devido a esta propriedade, denominada hierarquia, se pode construiu o Dendograma. Dendograma é um gráfico em forma de árvore no qual a escala vertical indica o nível de similaridade (ou dissimilaridade). No eixo horizontal tem-se os elementos e no vertical a altura correspondente ao nível em que os elementos foram considerados semelhantes. Quanto maior a altura maior a variabilidade, ou seja, maior a heterogeneidade dos grupos a serem unidos.

A outra técnica de construção dos clusters é chamada de Método hierárquivo Divisivo, que se inicia com todos os objetos em um grupo; formam-se subgrupos desagregando os grupos formados até terminar com cada objeto formando um grupo.

Existem ainda os Métodos não hierárquicos, também denominados de métodos de partição. Inicia-se com um número pré definido de grupos e a cada passo procura-se realocar os objetos de maneira a encontrar a melhor partição, isto é, a que minimiza a variância dentro do grupo e maximiza a variância entre grupos. Os métodos não hierárquicos não geram dendrograma, são aplicados a casos (observações) e não a variáveis e são mais indicados para amostras grandes.

2.4 Métodos de agrupamento

1. Método de ligação simples: também denominado de vizinho mais próximo ou single linkage. A similaridade entre dois grupos é definida pelos elementos mais parecidos. É a distância entre os vizinhos mais próximos ou entre os elementos mais parecidos de cada conglomerado.

\[ D_{I,II} = min \quad ({d_{ij}, i \in I, j \in II}) \]

2. Método de ligação completa ou do vizinho mais distante: A similaridade entre dois grupos é definida pelos elementos que são “menos semelhantes” entre si. Em cada passo do algoritmo são combinados os elementos que apresentarem o menor valor máximo de distância.

\[ D_{I,II} = max \quad ({d_{ij}, i \in I, j \in II}) \]

3. Método da média das distâncias (average linkage): Considera a distância entre dois grupos como a média das distâncias entre todos os pares de elementos dos dois grupos que estão sendo comparados.

\[ D_{i,j}=\frac{\sum_{i=1}^{n_I} \sum_{j=1}^{n_{II}}d_{ij}}{n_I n_{II}} \]

4. Método dos centróides: A distância entre dois grupos é definida como sendo a distância entre os vetores de médias (centróides), dos grupos comparados. É a distância Euclidiana quadrática entre os vetores de médias amostrais.

\[ D_{AB}=(\bar X_a - \bar X_b)'(\bar X_a - \bar X_b) \]

22. Método de Ward: No processo aglomerativo a medida que se agrupa o nível de similaridade diminui. Em cada passo do agrupamento ocorre diminuição de variabilidade entre grupos e aumento dentro dos grupos. O método de Ward se baseia na mudança de variação que ocorre de um estágio para outro. Este método tende a formar grupos com maior homogeneidade interna. O método considera o aumento na soma de quadrados dos erros como critério para juntar dois grupos. Para um dado grupo g, a soma de quadrados dos erros, \(SQE_g\), é a soma dos desvios ao quadrado de cada item para a média do grupo (centróide).

Em um estágio com G grupos, define-se \(SQE\) como a soma \(SQE=SQE_1+SQE_2+SQE_3+ \dots SQE_G\). Quando se junta dois grupos, SQE aumenta. Em cada etapa, a união de todos os pares de grupos é analisada. Junta-se os dois grupos que resultam no menor aumento de \(SQE\)

Segundo Mingoti (2005), “os métodos de ligação simples, completa e da média podem ser utilizados tanto para variáveis quantitativas, quanto qualitativas, ao contrário dos métodos do centróide e de Ward que são apropriados apenas para variáveis quantitativas, já que têm como base a comparação de vetores de médias”.

2.5 Análise de Conglomerados (Clusters) no R - Método Hieráquico

#Direcionado o R para o Diretorio a ser trabalhado
setwd('/Users/jricardofl/Dropbox/tempecon/dados_gescilene')

#Lendo os dados no R
library(foreign)
library(car)
library(tidyverse)
#library(Rcmdr)
library(corrplot)
library(MVar.pt)
library(cluster)
library(clValid) # Avaliação dos grupos
library(e1071) # Fuzzy K-médias
library(factoextra) # Vizualização de grupos
library(gridExtra) # Ferramentas gráficas
library(ggforce) # Ferramentas gráficas

#Entrada de Dados
dados <- read.dta("dados_gesc.dta")
str(dados)

## 'data.frame':    85 obs. of  16 variables:
##  $ x1 : num  60 56 56 61 47 65 50 60 78 68 ...
##  $ x2 : num  1 2 2 3 3 3 4 3 2 3 ...
##  $ x3 : num  140000 86100 47600 90000 205200 ...
##  $ x4 : num  15714 10500 9000 16667 17333 ...
##  $ x5 : num  17376 0 0 0 26064 ...
##  $ x6 : num  4502 1350 2400 2202 3376 ...
##  $ x7 : num  26400 21600 20400 18000 22800 21600 24000 13200 16800 21600 ...
##  $ x8 : num  48878 25800 26080 22902 55204 ...
##  $ x9 : num  2 0 0 0 3 0 2 0 2 0 ...
##  $ x10: num  0.111 0 0 0.222 0.111 ...
##  $ x11: num  0.143 0 0.143 0 0 ...
##  $ x12: num  20 0 1 5 20 20 10 15 20 10 ...
##  $ x13: num  2 0 1 2 2 2 2 1 1 1 ...
##  $ x14: num  1 0 1 0 0 0 0 0 0 0 ...
##  $ x15: num  0.26 0.12 0.09 0.06 0.06 ...
##  $ x16: num  1 0 1 0 1 0 1 1 1 1 ...
##  - attr(*, "datalabel")= chr ""
##  - attr(*, "time.stamp")= chr ""
##  - attr(*, "formats")= chr [1:16] "%9.0g" "%9.0g" "%9.0g" "%9.0g" ...
##  - attr(*, "types")= int [1:16] 254 254 254 254 254 254 254 254 254 254 ...
##  - attr(*, "val.labels")= chr [1:16] "" "" "" "" ...
##  - attr(*, "var.labels")= chr [1:16] "" "" "" "" ...
##  - attr(*, "version")= int -7

# Nomes das Variáveis utilizadas 
#x1 Idade
#x2 Escolaridade
#x3 Renda Bruta Total
#x4 Produtividade
#x5 Valor total anual  de mão de obra permanente
#x6 Valor Total com custo anual com insumos agrícolas nas ativ. Irrig. em 2014
#x7 Custo com energia eletrica e agua
#x8 Capital total empregado na ativ irrig em 2014
#x9 N total de empregados
#x10 Indice de introducao de inovacao (III)
#x11 Indice de Inovacoes realizadas (II)
#x12 Gastos para desenvolver as atividade de inovacao sobre faturamento de 2014
#x13 Quantidade de tecnologia agricola utilizada na atividade irrigada  
#x14 Dummy se a empresa realizou atividades de treinamento e capacitacao de recursos humanos entre 2010 a 2014
#x15 Indice de Fonte de Informacao
#x16 Dummy se empresa esteve envolvida em atividades cooperativas entre os anos de 2010 e 2014
#####################################################

#dados <- dados1 %>% select(-x14, -x16)
dados <- dados[,-c(14,16)]

#Se quiser padronizar os dados
dados.pad <-as.data.frame(scale(dados))

#Analise Fatorial
fatorial1 <- factanal(dados.pad, factors=4, rotation="none", na.action=na.omit)
fatorial2 <- factanal(dados.pad, factors=4, rotation="varimax", na.action=na.omit, scores = c("regression"))

#Analise de Cluster
# Função que obtém o WSS para o método hierárquico
# NÃO MUDAR!
get_wss <- function(d, cluster){
  d <- stats::as.dist(d)
  cn <- max(cluster)
  clusterf <- as.factor(cluster)
  clusterl <- levels(clusterf)
  cnn <- length(clusterl)
  
  if (cn != cnn) {
    warning("cluster renumbered because maximum != number of clusters")
    for (i in 1:cnn) cluster[clusterf == clusterl[i]] <- i
    cn <- cnn
  }
  cwn <- cn
  # Compute total within sum of square
  dmat <- as.matrix(d)
  within.cluster.ss <- 0
  for (i in 1:cn) {
    cluster.size <- sum(cluster == i)
    di <- as.dist(dmat[cluster == i, cluster == i])
    within.cluster.ss <- within.cluster.ss + sum(di^2)/cluster.size
  }
  within.cluster.ss
}

# Função que constrói o "gráfico do cotovelo" e aponta o K ótimo
# NÃO MUDAR!
elbow.plot <- function(x, kmax = 15, alg = "kmeans") {
  # alg = c("kmeans", "cmeans", "hclust")
  wss <- c()
  if (alg == "kmeans") {
    for (i in 1:kmax) {
      set.seed(13)
      tmp <- kmeans(x, i)
      # wss[i] <- get_wss(dist(x), tmp$cluster)
      wss[i] <- tmp$tot.withinss
    }
    tmp <- data.frame(k = 1:kmax, wss)
    max_k <- max(tmp$k)
    max_k_wss <- tmp$wss[which.max(tmp$k)]
    max_wss <- max(tmp$wss)
    max_wss_k <- tmp$k[which.max(tmp$wss)]
    max_df <- data.frame(x = c(max_wss_k, max_k), y = c(max_wss, max_k_wss))
    tmp_lm <- lm(max_df$y ~ max_df$x)
    d <- c()
    for(i in 1:kmax) {
      d <- c(d, abs(coef(tmp_lm)[2]*i - tmp$wss[i] + coef(tmp_lm)[1]) /
               sqrt(coef(tmp_lm)[2]^2 + 1^2))
    }
    tmp$d <- d
    ggplot(data = tmp, aes(k, wss)) +
      geom_line() +
      geom_segment(aes(x = k[1], y = wss[1],
                       xend = max(k), yend = wss[which.max(k)]),
                   linetype = "dashed") +
      geom_point(aes(size = (d == max(d)), color = (d == max(d))),
                 show.legend = FALSE) +
      scale_size_manual(values = c(2,5)) +
      scale_color_manual(values = c("black", "red")) +
      labs(x = "Number of clusters",
           y = "Total within-cluster sum of squares",
           title = "Elbow plot for the K-means method") +
      theme_bw()
  }
  else if (alg == "cmeans") {
    for (i in 1:kmax) {
      if (i == 1) {
        wss[i] <- get_wss(dist(x), rep(1, nrow(x)))
      }
      else {
        set.seed(13)
        tmp <- cmeans(x, i)
        wss[i] <- get_wss(dist(x), tmp$cluster)
        # wss[i] <- tmp$sumsqrs$tot.within.ss
      }
    }
    tmp <- data.frame(k = 1:kmax, wss)
    max_k <- max(tmp$k)
    max_k_wss <- tmp$wss[which.max(tmp$k)]
    max_wss <- max(tmp$wss)
    max_wss_k <- tmp$k[which.max(tmp$wss)]
    max_df <- data.frame(x = c(max_wss_k, max_k), y = c(max_wss, max_k_wss))
    tmp_lm <- lm(max_df$y ~ max_df$x)
    d <- c()
    for(i in 1:kmax) {
      d <- c(d, abs(coef(tmp_lm)[2]*i - tmp$wss[i] + coef(tmp_lm)[1]) /
               sqrt(coef(tmp_lm)[2]^2 + 1^2))
    }
    tmp$d <- d
    ggplot(data = tmp, aes(k, wss)) +
      geom_line() +
      geom_segment(aes(x = k[1], y = wss[1],
                       xend = max(k), yend = wss[which.max(k)]),
                   linetype = "dashed") +
      geom_point(aes(size = (d == max(d)), color = (d == max(d))),
                 show.legend = FALSE) +
      scale_size_manual(values = c(2,5)) +
      scale_color_manual(values = c("black", "red")) +
      labs(x = "Number of clusters",
           y = "Total within-cluster sum of squares",
           title = "Elbow plot for the fuzzy K-means method") +
      theme_bw()
  }
  else if (alg == "hclust") {
    for (i in 1:kmax) {
      set.seed(13)
      tmp <- hcut(x, i)
      wss[i] <- get_wss(dist(x), tmp$cluster)
    }
    tmp <- data.frame(k = 1:kmax, wss)
    max_k <- max(tmp$k)
    max_k_wss <- tmp$wss[which.max(tmp$k)]
    max_wss <- max(tmp$wss)
    max_wss_k <- tmp$k[which.max(tmp$wss)]
    max_df <- data.frame(x = c(max_wss_k, max_k), y = c(max_wss, max_k_wss))
    tmp_lm <- lm(max_df$y ~ max_df$x)
    d <- c()
    for(i in 1:kmax) {
      d <- c(d, abs(coef(tmp_lm)[2]*i - tmp$wss[i] + coef(tmp_lm)[1]) /
               sqrt(coef(tmp_lm)[2]^2 + 1^2))
    }
    tmp$d <- d
    ggplot(data = tmp, aes(k, wss)) +
      geom_line() +
      geom_segment(aes(x = k[1], y = wss[1],
                       xend = max(k), yend = wss[which.max(k)]),
                   linetype = "dashed") +
      geom_point(aes(size = (d == max(d)), color = (d == max(d))),
                 show.legend = FALSE) +
      scale_size_manual(values = c(2,5)) +
      scale_color_manual(values = c("black", "red")) +
      labs(x = "Number of clusters",
           y = "Total within-cluster sum of squares",
           title = "Elbow plot for the hierarchical method") +
      theme_bw()
  }
}

# Função vizualização dos grupos
# NÃO MUDAR!
cluster_viz <- function(data, clusters, 
                        axes = c(1, 2), geom = c("point", "text"), repel = TRUE, 
                        show.clust.cent = TRUE, ellipse = TRUE, ellipse.type = "convex", 
                        ellipse.level = 0.95, ellipse.alpha = 0.2, shape = NULL, 
                        pointsize = 1.5, labelsize = 12, main = "Cluster plot",
                        ggtheme = theme_bw()) {
  require(factoextra)
  data <- scale(data)
  pca <- stats::prcomp(data, scale = FALSE, center = FALSE)
  ind <- facto_summarize(pca, element = "ind", result = "coord", axes = axes)
  eig <- get_eigenvalue(pca)[axes, 2]
  xlab <- paste0("Dim", axes[1], " (", round(eig[1], 1), "%)")
  ylab <- paste0("Dim", axes[2], " (", round(eig[2], 1), "%)")
  colnames(ind)[2:3] <- c("x", "y")
  label_coord <- ind
  lab <- NULL
  if ("text" %in% geom) 
    lab <- "name"
  if (is.null(shape)) 
    shape <- "cluster"
  plot.data <- cbind.data.frame(ind, cluster = clusters, stringsAsFactors = TRUE)
  label_coord <- cbind.data.frame(label_coord, cluster = clusters, stringsAsFactors = TRUE)
  p <- ggpubr::ggscatter(plot.data, "x", "y", color = "cluster", 
                         shape = shape, size = pointsize, point = "point" %in% geom,
                         label = lab, font.label = labelsize, repel = repel, 
                         mean.point = show.clust.cent, ellipse = ellipse, ellipse.type = ellipse.type, 
                         ellipse.alpha = ellipse.alpha, ellipse.level = ellipse.level, 
                         main = main, xlab = xlab, ylab = ylab, ggtheme = ggtheme)
  p
}

# Encontrar número ótimo de grupos para o método hierárquico
set.seed(123456789)
p3 <- elbow.plot(dados.pad, alg = "hclust")
p3

# K "ótimo" é igual a 4. Mudar de acordo com o seu resultado

d <-dist(dados.pad, method = "euclidean")
h.fit <-hclust(d, method = "ward.D")
plot(h.fit, main='Dendograma - Método Hierárquico', xlab='Cluster das Observações - Distância Euclidiana e Método de Ward', ylab='Altura')

groups_a <- cutree(h.fit, k=3)
groups_a 

##  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 
##  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1 
## 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 
##  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  2  1  2 
## 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 
##  2  2  1  2  2  2  2  2  2  1  2  2  2  2  2  2  2  3  2  3  2  2  3  2  2  3 
## 79 80 81 82 83 84 85 
##  3  2  2  2  2  2  2

groups_b <- cutree(h.fit, k=4)
groups_b 

##  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 
##  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1 
## 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 
##  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  2  1  3 
## 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 
##  3  3  1  2  2  2  2  2  3  1  2  2  3  3  2  2  2  4  2  4  2  2  4  3  2  4 
## 79 80 81 82 83 84 85 
##  4  2  2  2  2  2  2

groups_c <- cutree(h.fit, k=5)
dendogram3 <- rect.hclust(h.fit, k=5, border="red")

dendogram3

## [[1]]
##  2  3  4  6  9 10 11 13 17 18 28 29 31 33 37 42 43 48 49 55 62 
##  2  3  4  6  9 10 11 13 17 18 28 29 31 33 37 42 43 48 49 55 62 
## 
## [[2]]
##  1  5  7  8 12 14 15 16 19 20 21 22 23 24 25 26 27 30 32 34 35 36 38 39 40 41 
##  1  5  7  8 12 14 15 16 19 20 21 22 23 24 25 26 27 30 32 34 35 36 38 39 40 41 
## 44 45 46 47 51 
## 44 45 46 47 51 
## 
## [[3]]
## 70 72 75 78 79 
## 70 72 75 78 79 
## 
## [[4]]
## 52 53 54 61 65 66 76 
## 52 53 54 61 65 66 76 
## 
## [[5]]
## 50 56 57 58 59 60 63 64 67 68 69 71 73 74 77 80 81 82 83 84 85 
## 50 56 57 58 59 60 63 64 67 68 69 71 73 74 77 80 81 82 83 84 85

summary(fatorial2$scores)

##     Factor1           Factor2            Factor3           Factor4        
##  Min.   :-1.1249   Min.   :-1.11082   Min.   :-3.1932   Min.   :-2.11638  
##  1st Qu.:-0.6559   1st Qu.:-0.21990   1st Qu.:-0.2904   1st Qu.:-0.18303  
##  Median :-0.3633   Median :-0.11821   Median :-0.1839   Median :-0.14019  
##  Mean   : 0.0000   Mean   : 0.00000   Mean   : 0.0000   Mean   : 0.00000  
##  3rd Qu.: 0.3353   3rd Qu.:-0.04689   3rd Qu.: 0.0107   3rd Qu.:-0.06075  
##  Max.   : 3.3937   Max.   : 7.06890   Max.   : 6.5278   Max.   : 7.41904

attach(dados)
by(dados, groups_c, summary)

## groups_c: 1
##        x1             x2              x3               x4       
##  Min.   :29.0   Min.   :1.000   Min.   : 33000   Min.   : 6667  
##  1st Qu.:49.0   1st Qu.:2.000   1st Qu.: 87150   1st Qu.:13345  
##  Median :52.0   Median :3.000   Median :127400   Median :16667  
##  Mean   :53.1   Mean   :2.903   Mean   :147082   Mean   :16839  
##  3rd Qu.:60.0   3rd Qu.:3.000   3rd Qu.:166400   3rd Qu.:21027  
##  Max.   :70.0   Max.   :5.000   Max.   :461700   Max.   :30000  
##        x5              x6              x7              x8        
##  Min.   :    0   Min.   : 1071   Min.   :12000   Min.   : 16660  
##  1st Qu.:    0   1st Qu.: 1401   1st Qu.:18000   1st Qu.: 30268  
##  Median :17376   Median : 1900   Median :19200   Median : 43008  
##  Mean   :18577   Mean   : 3339   Mean   :20923   Mean   : 44942  
##  3rd Qu.:26064   3rd Qu.: 3378   3rd Qu.:24000   3rd Qu.: 55152  
##  Max.   :69504   Max.   :27480   Max.   :38400   Max.   :116484  
##        x9             x10              x11               x12      
##  Min.   :0.000   Min.   :0.0000   Min.   :0.00000   Min.   :10.0  
##  1st Qu.:0.000   1st Qu.:0.1111   1st Qu.:0.00000   1st Qu.:20.0  
##  Median :2.000   Median :0.1111   Median :0.00000   Median :20.0  
##  Mean   :2.065   Mean   :0.1254   Mean   :0.03226   Mean   :22.9  
##  3rd Qu.:3.000   3rd Qu.:0.1111   3rd Qu.:0.00000   3rd Qu.:25.0  
##  Max.   :8.000   Max.   :0.2222   Max.   :0.14286   Max.   :50.0  
##       x13             x15      
##  Min.   :0.000   Min.   :0.03  
##  1st Qu.:1.000   1st Qu.:0.06  
##  Median :1.000   Median :0.12  
##  Mean   :1.323   Mean   :0.12  
##  3rd Qu.:2.000   3rd Qu.:0.17  
##  Max.   :2.000   Max.   :0.30  
## ------------------------------------------------------------ 
## groups_c: 2
##        x1              x2              x3               x4       
##  Min.   :50.00   Min.   :1.000   Min.   : 21250   Min.   : 3333  
##  1st Qu.:61.00   1st Qu.:2.000   1st Qu.: 49200   1st Qu.:13333  
##  Median :65.00   Median :2.000   Median : 86100   Median :16429  
##  Mean   :65.95   Mean   :2.095   Mean   : 91014   Mean   :17585  
##  3rd Qu.:71.00   3rd Qu.:2.000   3rd Qu.:126000   3rd Qu.:21000  
##  Max.   :78.00   Max.   :3.000   Max.   :160402   Max.   :31000  
##        x5              x6              x7              x8       
##  Min.   :    0   Min.   :  880   Min.   : 9240   Min.   :12870  
##  1st Qu.:    0   1st Qu.: 1650   1st Qu.:16800   1st Qu.:23410  
##  Median :    0   Median : 2003   Median :20400   Median :26080  
##  Mean   : 8834   Mean   : 3469   Mean   :19069   Mean   :33256  
##  3rd Qu.:17376   3rd Qu.: 3000   3rd Qu.:21600   3rd Qu.:38000  
##  Max.   :34752   Max.   :25300   Max.   :26400   Max.   :62740  
##        x9              x10              x11               x12        
##  Min.   :0.0000   Min.   :0.0000   Min.   :0.00000   Min.   : 0.000  
##  1st Qu.:0.0000   1st Qu.:0.1111   1st Qu.:0.00000   1st Qu.: 5.000  
##  Median :0.0000   Median :0.1111   Median :0.00000   Median :10.000  
##  Mean   :0.9524   Mean   :0.1058   Mean   :0.03401   Mean   : 8.619  
##  3rd Qu.:2.0000   3rd Qu.:0.1111   3rd Qu.:0.00000   3rd Qu.:10.000  
##  Max.   :4.0000   Max.   :0.2222   Max.   :0.14286   Max.   :20.000  
##       x13             x15       
##  Min.   :0.000   Min.   :0.060  
##  1st Qu.:1.000   1st Qu.:0.120  
##  Median :1.000   Median :0.120  
##  Mean   :1.286   Mean   :0.149  
##  3rd Qu.:2.000   3rd Qu.:0.160  
##  Max.   :2.000   Max.   :0.460  
## ------------------------------------------------------------ 
## groups_c: 3
##        x1              x2              x3                x4       
##  Min.   :26.00   Min.   :1.000   Min.   : 105000   Min.   :11284  
##  1st Qu.:35.00   1st Qu.:5.000   1st Qu.: 208000   1st Qu.:22381  
##  Median :42.00   Median :5.000   Median : 403580   Median :30000  
##  Mean   :43.24   Mean   :4.952   Mean   : 567728   Mean   :30888  
##  3rd Qu.:51.00   3rd Qu.:6.000   3rd Qu.: 720000   3rd Qu.:34714  
##  Max.   :65.00   Max.   :8.000   Max.   :2380000   Max.   :62963  
##        x5               x6               x7              x8        
##  Min.   : 10080   Min.   :  1350   Min.   : 3720   Min.   : 51280  
##  1st Qu.: 42000   1st Qu.: 55000   1st Qu.:12360   1st Qu.:148600  
##  Median : 83040   Median :145000   Median :15840   Median :309824  
##  Mean   :117334   Mean   :146283   Mean   :20742   Mean   :306951  
##  3rd Qu.:167202   3rd Qu.:220000   3rd Qu.:25800   3rd Qu.:420600  
##  Max.   :312000   Max.   :350000   Max.   :51600   Max.   :608400  
##        x9              x10              x11              x12        
##  Min.   : 0.000   Min.   :0.0000   Min.   :0.0000   Min.   : 0.000  
##  1st Qu.: 4.000   1st Qu.:0.2222   1st Qu.:0.0000   1st Qu.: 0.000  
##  Median : 8.000   Median :0.3333   Median :0.1429   Median : 5.000  
##  Mean   : 9.238   Mean   :0.2910   Mean   :0.1463   Mean   : 5.048  
##  3rd Qu.:13.000   3rd Qu.:0.4444   3rd Qu.:0.2857   3rd Qu.:10.000  
##  Max.   :25.000   Max.   :0.5556   Max.   :0.4286   Max.   :10.000  
##       x13             x15        
##  Min.   :1.000   Min.   :0.0000  
##  1st Qu.:3.000   1st Qu.:0.4700  
##  Median :4.000   Median :0.5900  
##  Mean   :3.571   Mean   :0.5671  
##  3rd Qu.:4.000   3rd Qu.:0.7100  
##  Max.   :6.000   Max.   :0.8900  
## ------------------------------------------------------------ 
## groups_c: 4
##        x1              x2              x3                x4       
##  Min.   :39.00   Min.   :3.000   Min.   : 176000   Min.   : 8000  
##  1st Qu.:49.50   1st Qu.:6.000   1st Qu.: 714000   1st Qu.:13583  
##  Median :52.00   Median :7.000   Median :1015000   Median :19444  
##  Mean   :52.43   Mean   :6.143   Mean   :1017429   Mean   :19537  
##  3rd Qu.:57.00   3rd Qu.:7.000   3rd Qu.:1446500   3rd Qu.:20800  
##  Max.   :63.00   Max.   :7.000   Max.   :1610000   Max.   :40545  
##        x5               x6               x7              x8        
##  Min.   : 24000   Min.   : 49000   Min.   :15240   Min.   :101280  
##  1st Qu.:125544   1st Qu.:134150   1st Qu.:28380   1st Qu.:472710  
##  Median :242448   Median :235000   Median :37200   Median :605008  
##  Mean   :223483   Mean   :244614   Mean   :44286   Mean   :588103  
##  3rd Qu.:279696   3rd Qu.:330000   3rd Qu.:59400   3rd Qu.:739880  
##  Max.   :487452   Max.   :500000   Max.   :82000   Max.   :985252  
##        x9             x10              x11              x12       
##  Min.   : 0.00   Min.   :0.5556   Min.   :0.4286   Min.   : 1.00  
##  1st Qu.: 8.50   1st Qu.:0.5556   1st Qu.:0.5000   1st Qu.: 6.50  
##  Median :22.00   Median :0.5556   Median :0.6429   Median :10.00  
##  Mean   :20.14   Mean   :0.6508   Mean   :0.5918   Mean   :11.71  
##  3rd Qu.:28.50   3rd Qu.:0.7222   3rd Qu.:0.6786   3rd Qu.:15.00  
##  Max.   :45.00   Max.   :0.8889   Max.   :0.7143   Max.   :28.00  
##       x13              x15        
##  Min.   : 5.000   Min.   :0.6700  
##  1st Qu.: 6.000   1st Qu.:0.7800  
##  Median : 6.000   Median :0.8800  
##  Mean   : 7.286   Mean   :0.8629  
##  3rd Qu.: 7.500   3rd Qu.:0.9650  
##  Max.   :13.000   Max.   :1.0000  
## ------------------------------------------------------------ 
## groups_c: 5
##        x1             x2            x3                x4       
##  Min.   :27.0   Min.   :7.0   Min.   :4290000   Min.   :18033  
##  1st Qu.:29.0   1st Qu.:7.0   1st Qu.:4375000   1st Qu.:25000  
##  Median :29.0   Median :7.0   Median :6160000   Median :34310  
##  Mean   :35.4   Mean   :7.4   Mean   :5616500   Mean   :29814  
##  3rd Qu.:40.0   3rd Qu.:8.0   3rd Qu.:6237500   3rd Qu.:35227  
##  Max.   :52.0   Max.   :8.0   Max.   :7020000   Max.   :36500  
##        x5                x6                x7               x8         
##  Min.   : 391680   Min.   : 553000   Min.   :     0   Min.   :1829004  
##  1st Qu.: 456840   1st Qu.:1200000   1st Qu.: 60000   1st Qu.:1855840  
##  Median : 920004   Median :1300000   Median : 60000   Median :2747440  
##  Mean   : 809705   Mean   :1730600   Mean   :105600   Mean   :2991657  
##  3rd Qu.:1080000   3rd Qu.:1900000   3rd Qu.:144000   3rd Qu.:3376000  
##  Max.   :1200000   Max.   :3700000   Max.   :264000   Max.   :5150000  
##        x9            x10              x11              x12          x13      
##  Min.   : 0.0   Min.   :0.4444   Min.   :0.2143   Min.   : 0   Min.   : 3.0  
##  1st Qu.:30.0   1st Qu.:0.7778   1st Qu.:0.4286   1st Qu.: 5   1st Qu.: 6.0  
##  Median :44.0   Median :0.7778   Median :0.7143   Median : 5   Median : 8.0  
##  Mean   :50.8   Mean   :0.7556   Mean   :0.6143   Mean   : 8   Mean   : 7.4  
##  3rd Qu.:90.0   3rd Qu.:0.7778   3rd Qu.:0.8571   3rd Qu.:10   3rd Qu.: 8.0  
##  Max.   :90.0   Max.   :1.0000   Max.   :0.8571   Max.   :20   Max.   :12.0  
##       x15       
##  Min.   :0.460  
##  1st Qu.:0.760  
##  Median :0.810  
##  Mean   :0.762  
##  3rd Qu.:0.860  
##  Max.   :0.920

# Agrupamento usando o método hierárquico - 2 forma
nclust = 5 # Número de grupos de acordo com os gráficos
fit.hclust <- hcut(dados.pad, 5) # Ajustar K de acordo com o gráfico
fit.hclust

## 
## Call:
## stats::hclust(d = x, method = hc_method)
## 
## Cluster method   : ward.D2 
## Distance         : euclidean 
## Number of objects: 85

g_cluster <- cluster_viz(dados.pad, as.factor(fit.hclust$cluster), geom = "point", main = "Gráfico Cluster para o Método Hierárquico")
plot(g_cluster)

2.6 Métodos Não-Hierárquicos

São métodos de partição de “n” objetos em “k” grupos. Os “k” grupos são pré-determinados. A cada passo verifica-se se os objetos estão alocados da “melhor” maneira. Função objetivo: procura minimizar a variância dentro do grupo e maximizar a variância entre os grupos.

Não fornece dendograma. Isto porque se em algum passo do algoritmo dois elementos tiverem sido colocados num mesmo cluster, não necessariamente “estarão juntos” na partição final. Bom quando se tem um grande número de observações. O problema principal aqui não é o número de grupos e sim a melhor forma de alocar os elementos em k grupos. O método mais usado é o das k-médias (k-means).

Resumidamente, o método das k-means é composto de quatro passos. No primeiro são selecionados k centróides para inicializar o processo de partição. Cada elemento do conjunto de dados é, então, comparado com cada centróide inicial, através de uma medida de distância que, em geral, é a Euclidiana. O elemento é alocado ao grupo cuja distância é a menor. No terceiro passo, recalcula-se os valores dos centróides para cada novo grupo formado. Com estes novos centróides, repete-se o passo 2. Isto deve ser repetido até que nenhuma realocação de elementos seja necessária.

cluster_kmeans <- elbow.plot(dados.pad)
cluster_kmeans

kmeans1 <- kmeans(dados.pad, centers=5)
#Centróides
kmeans1$centers

##            x1         x2           x3         x4         x5          x6
## 1 -1.34515006  1.8072111  3.701528552  0.8027013  3.2778778  3.23549118
## 2 -0.02424868  1.1862676  0.287208385 -0.1871013  0.5550533  0.17399627
## 3 -0.02004540 -0.3831231 -0.355117495 -0.3418032 -0.3879402 -0.32287084
## 4 -0.79009568  0.6492143  0.007068955  1.1307878  0.1277874  0.01285434
## 5  0.83455884 -0.6890467 -0.397959050 -0.5056032 -0.4355642 -0.31782898
##           x7          x8         x9        x10        x11        x12        x13
## 1  2.4501995  3.49247328  2.7059130  2.2639164  2.1490201 -0.5026098  1.8761180
## 2  0.5309120  0.34475417  0.7714212  1.8022454  2.0470493 -0.1411494  1.8303262
## 3 -0.1966927 -0.36414512 -0.3560014 -0.5085221 -0.5070145  1.0175370 -0.5500592
## 4 -0.1810064  0.02647696  0.1454207  0.2944307  0.1338101 -0.7891520  0.3802517
## 5 -0.2694298 -0.37661625 -0.4316555 -0.5572274 -0.4915234 -0.4539516 -0.5032957
##          x15
## 1  1.4079157
## 2  1.7417060
## 3 -0.6848525
## 4  0.9214138
## 5 -0.6137018

#Plot dos Clusters
clusplot(dados.pad, kmeans1$cluster, main="2D Representação da solução dos Clusters",
         color=TRUE, shade=TRUE, labels=2, lines=0)

# Agrupamento usando o K-means - 2 forma
fit.kmeans <- kmeans(dados.pad, 5) # Ajustar K de acordo com o gráfico
fit.kmeans

## K-means clustering with 5 clusters of sizes 18, 5, 29, 26, 7
## 
## Cluster means:
##            x1         x2           x3         x4         x5          x6
## 1 -0.79009568  0.6492143  0.007068955  1.1307878  0.1277874  0.01285434
## 2 -1.34515006  1.8072111  3.701528552  0.8027013  3.2778778  3.23549118
## 3 -0.02004540 -0.3831231 -0.355117495 -0.3418032 -0.3879402 -0.32287084
## 4  0.83455884 -0.6890467 -0.397959050 -0.5056032 -0.4355642 -0.31782898
## 5 -0.02424868  1.1862676  0.287208385 -0.1871013  0.5550533  0.17399627
##           x7          x8         x9        x10        x11        x12        x13
## 1 -0.1810064  0.02647696  0.1454207  0.2944307  0.1338101 -0.7891520  0.3802517
## 2  2.4501995  3.49247328  2.7059130  2.2639164  2.1490201 -0.5026098  1.8761180
## 3 -0.1966927 -0.36414512 -0.3560014 -0.5085221 -0.5070145  1.0175370 -0.5500592
## 4 -0.2694298 -0.37661625 -0.4316555 -0.5572274 -0.4915234 -0.4539516 -0.5032957
## 5  0.5309120  0.34475417  0.7714212  1.8022454  2.0470493 -0.1411494  1.8303262
##          x15
## 1  0.9214138
## 2  1.4079157
## 3 -0.6848525
## 4 -0.6137018
## 5  1.7417060
## 
## Clustering vector:
##  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 
##  3  4  4  4  3  3  4  4  4  4  4  3  4  3  3  3  4  4  3  3  3  3  3  3  3  3 
## 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 
##  3  4  4  4  4  3  4  3  3  3  4  3  3  3  4  4  4  3  3  3  3  4  4  3  3  5 
## 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 
##  5  5  4  1  4  1  1  1  5  4  1  1  5  5  1  1  1  2  1  2  1  1  2  5  1  2 
## 79 80 81 82 83 84 85 
##  2  1  4  1  1  1  1 
## 
## Within cluster sum of squares by cluster:
## [1]  86.28616 143.25082  64.37370  48.86236  43.95670
##  (between_SS / total_SS =  67.1 %)
## 
## Available components:
## 
## [1] "cluster"      "centers"      "totss"        "withinss"     "tot.withinss"
## [6] "betweenss"    "size"         "iter"         "ifault"

g_kmeans <- cluster_viz(dados.pad, as.factor(fit.kmeans$cluster), geom = "point",
                  main = "Gráfico Cluster para o Método K-médias")
plot(g_kmeans)

3 Análise Discriminante

A Análise Discriminante Linear (LDA) é uma técnica de Análise Multivariada usada com o objetivo de discriminar (separar, diferenciar) populações e/ou classificar observações em populações pré definidas. Para a sua aplicação, é necessário que os grupos para os quais cada elemento amostral pode ser classificado sejam pré-definidos, ou seja, conhecidos a priori considerando-se suas características gerais. Este conhecimento permite a elaboração de uma função matemática chamada de regra de classificação ou discriminação, que é utilizada para classificar novos elementos amostrais nos grupos já existentes. A definição dos grupos é feita de acordo com o problema de pesquisa e os objetivos do estudo.

A LDA, originalmente desenvolvida por R A Fisher (1936) para classificar indivíduos em um dos dois grupos claramente definidos. Posteriormente, foi expandido para classificar os assuntos em mais de dois grupos. A LDA é uma técnica para reduzir o número de dimensões (ou seja, variáveis) em um conjunto de dados, mantendo o máximo de informações possível. Basicamente, ajuda a encontrar a combinação linear das variáveis originais que fornecem a melhor separação possível entre os grupos.

O objetivo é construir uma regra de classificação, baseada na teoria das probabilidades, que minimize o número de classificações incorretas, ou seja, o erro de dizer que um elemento pertence à população A quando na verdade ele pertence à população B. Uma informação importante é saber quais variáveis mais ajudam a explicar a diferença entre os grupos. A idéia básica é saber se as respectivas médias dos grupos provêm de distintas populações.

A Análise Discriminante tem, assim, dois objetivos básicos:

1. Discriminar – consiste em encontrar funções capazes de explicar as diferenças entre os grupos com base em várias variáveis; estas funções captam o máximo de separação entre os grupos.

2. Classificar – consiste no uso das funções para alocar novas observações em um dos grupos pré-estabelecidos. Assim, pode-se dizer que existem dois tipos de Análise Discriminante: Descritiva e Preditiva.

As principais aplicações ou exemplos são:

1. Prevendo o sucesso ou fracasso de novos produtos;

2. Aceitar ou rejeitar a admissão de um candidato;

3. Prevendo a categoria de risco de crédito para uma pessoa;

4. Classificando os pacientes em diferentes categorias.

Para aplicação da LDA deve-se dispor de uma amostra de cada uma das g populações e dados de p variáveis observadas para cada objeto das amostras que podem ser de tamanhos diferentes. Ao contrário da Análise de Cluster em que se procura formar grupos homogêneos na amostra, na LDA parte-se de um número pré-definido de grupos considerados amostras de diferentes populações.

A base de dados utilizada para o exemplo é o da flor Iris ou conjunto de dados Iris de Fisher, um conjunto de dados multivariados introduzido pelo estatístico e biólogo britânico Ronald Fisher em seu artigo de 1936, “O uso de múltiplas medições em problemas taxonômicos, como um exemplo de análise discriminante linear”. Às vezes, é chamado de conjunto de dados da íris de Anderson porque Edgar Anderson coletou os dados para quantificar a variação morfológica das flores da íris de três espécies relacionadas. Duas das três espécies foram coletadas na Península de Gaspé, “todas do mesmo campo, colhidas no mesmo dia e medidas ao mesmo tempo pela mesma pessoa com a mesma aparelho”. O conjunto de dados consiste em 50 amostras de cada uma das três espécies de Iris ( Iris setosa, Iris virginica e Iris versicolor). Quatro variáveis foram medidas em cada amostra: o comprimento e a largura das sépalas e pétalas, em centímetros. Com base na combinação dessas quatro características, Fisher desenvolveu um modelo discriminante linear para distinguir as espécies umas das outras (Fonte: https://pt.wikipedia.org/wiki/Conjunto_de_dados_flor_Iris).

3.1 Análise Discriminante Linear no R

#Pacotes utilizados
library(klaR)
library(psych)
library(MASS)
library(ggord)
library(devtools)

# Obtendo os dados
data("iris")
str(iris)

## 'data.frame':    150 obs. of  5 variables:
##  $ Sepal.Length: num  5.1 4.9 4.7 4.6 5 5.4 4.6 5 4.4 4.9 ...
##  $ Sepal.Width : num  3.5 3 3.2 3.1 3.6 3.9 3.4 3.4 2.9 3.1 ...
##  $ Petal.Length: num  1.4 1.4 1.3 1.5 1.4 1.7 1.4 1.5 1.4 1.5 ...
##  $ Petal.Width : num  0.2 0.2 0.2 0.2 0.2 0.4 0.3 0.2 0.2 0.1 ...
##  $ Species     : Factor w/ 3 levels "setosa","versicolor",..: 1 1 1 1 1 1 1 1 1 1 ...

#Criando um gráfico de dispersão para as quatro primeiras variáveis numéricas.

pairs.panels(iris[1:4],
             gap = 0,
             bg = c("red", "green", "blue")[iris$Species],
             pch = 21)

Se pode ver no gráfico, diagrama de dispersão, histograma e valores de correlação. Em seguida, se quer criar os melhores grupos de separação com base nessas espécies. Vai ser criado um conjunto de dados de treinamento (training) e um conjunto de dados de teste (test) para fins de previsão e teste. 60% do conjunto de dados usado para fins de treinamento e 40% usados para fins de teste.

#Separação dos grupos
set.seed(123456789)
ind <- sample(2, nrow(iris),
              replace = TRUE,
              prob = c(0.6, 0.4))
training <- iris[ind==1,]
testing <- iris[ind==2,]

# Análise discriminante linear
linear <- lda(Species~., data=training)
linear

## Call:
## lda(Species ~ ., data = training)
## 
## Prior probabilities of groups:
##     setosa versicolor  virginica 
##  0.3043478  0.3478261  0.3478261 
## 
## Group means:
##            Sepal.Length Sepal.Width Petal.Length Petal.Width
## setosa         5.035714    3.460714     1.475000   0.2357143
## versicolor     5.896875    2.731250     4.215625   1.3125000
## virginica      6.553125    3.034375     5.509375   2.0281250
## 
## Coefficients of linear discriminants:
##                     LD1        LD2
## Sepal.Length  0.5851444  0.4528399
## Sepal.Width   2.0168111 -2.5347510
## Petal.Length -1.8643109  0.6817960
## Petal.Width  -3.2514946 -2.7190523
## 
## Proportion of trace:
##    LD1    LD2 
## 0.9855 0.0145

attributes(linear)

## $names
##  [1] "prior"   "counts"  "means"   "scaling" "lev"     "svd"     "N"      
##  [8] "call"    "terms"   "xlevels"
## 
## $class
## [1] "lda"

# Histograma empilhado para valores de funções discriminantes.
p <- predict(linear, training)
ldahist(data = p$x[,1], g = training$Species)

Com base no conjunto de dados de treinamento, 30,4% pertencem ao grupo setosa, 34,7% pertencem aos grupos versicolor e 34,7% pertencem aos grupos virginica. A primeira função discriminante é uma combinação linear das quatro variáveis. A porcentagem de separações alcançada pela primeira função discriminante é 98,55% e a segunda é 1,45%.

Esses histogramas são baseados na primeira função discriminante. É evidente que não há sobreposições entre primeira e segunda e primeira e terceira espécies. Contudo, alguma sobreposição é observada entre a segunda e a terceira espécie.

ldahist(data = p$x[,2], g = training$Species)

Os histogramas baseados na segunda função discriminante são sobrepostos, não bons.

ggord(linear, training$Species, ylim = c(-10, 10))

Os gráficos Biplot acima são baseados nas duas funções discriminantes. Setosa se separou muito claramente e alguma sobreposição foi observada entre Versicolor e Virginica.

Com base nas setas, a largura e o comprimento da sépala explicaram mais para Setosa, a largura e o comprimento da pétala explicaram mais para Versicolor e Virginica.

Finalmente, é possível ver as classificações corretas e classificações de erros.

p1 <- predict(linear, training)$class
tab <- table(Predicted = p1, Actual = training$Species)
tab

##             Actual
## Predicted    setosa versicolor virginica
##   setosa         28          0         0
##   versicolor      0         30         1
##   virginica       0          2        31

sum(diag(tab))/sum(tab)

## [1] 0.9673913

p2 <- predict(linear, testing)$class
tab1 <- table(Predicted = p2, Actual = testing$Species)
tab1

##             Actual
## Predicted    setosa versicolor virginica
##   setosa         22          0         0
##   versicolor      0         18         0
##   virginica       0          0        18

sum(diag(tab1))/sum(tab1)

## [1] 1

No conjunto de dados de treinamento, a classificação correta total é 30+30+29=89. A acurácia do modelo é aproximadamente 100%.

No conjunto de dados de teste, a classificação correta total é 20+19+20=59. A acurácia do modelo é aproximadamente 96,72%.

Assim, o Histograma e o Biplot fornecem informações úteis para interpretações e, se não houver uma grande diferença nas matrizes de covariância do grupo, a análise discriminante linear funcionará muito bem. A Análise Discriminante Linear não é útil para resolver problemas não lineares.